# -*- coding: utf-8 -*-
"""
Creating a normalised model from the original model

@author: obiri
"""
import mpctools as mpc
import numpy as np
import casadi as cs
import time

##16 states, 21 parameters. 
#nb: some of the remaining parameters are fixed measurements, one is state dependent

lb = 0.5 # lower bound of normalized model: lb*xs
ub = 1.5 # upper bound of normalized model: ub*xs
width = ub-lb

# xs = np.array([3.27367434e+10, 8.28831113e+10, 2.71741538e+03, 5.78604334e+00,
#         6.37656737e+01, 4.83137928e+00, 6.23866477e+01, 3.00039503e+03,
#         2.61960613e+09, 6.63233551e+09, 2.17448574e+03, 4.63001301e+00,
#         5.10255242e+01, 3.86608729e+00, 4.99220220e+01, 2.97694060e+03,
#         3.70349577e+01, 1.50000000e+01, 4.99220218e+01])

xs = np.array([3.27367434e+10, 8.28831113e+10, 2.71741538e+03, 5.78604334e+00,
        6.37656737e+01, 4.83137928e+00, 6.23866477e+01, 
        2.61960613e+09, 6.63233551e+09, 2.17448574e+03, 4.63001301e+00,
        5.10255242e+01, 3.86608729e+00, 4.99220220e+01, 
        3.70349577e+01, 4.99220218e+01])

def reactor_tank(x, u, w):   
    
    """
    An integrated model including upstream and buffer tank only.
    :param x: State vector (19,1)
    :param u: Input vector (8,1)
    :return: ODEs (19,1)
    """
    dxdt = cs.SX.zeros(16)
    
    # # States:
    # Xv1 = x[0]   # concentration of viable cells in reactor(cell/L)
    # Xt1 = x[1]   # total concentration of cells in reactor(cell/L)
    # GLC1 = x[2]  # glucose concentration in reactor(mM)
    # GLN1 = x[3]  # glutamine concentration in reactor(mM)
    # LAC1 = x[4]  # lactate concentration in reactor(mM)
    # AMM1 = x[5]  # ammonia concentration in reactor(mM)
    # mAb1 = x[6]  # mAb concentration in reactor(mg/L)
   ## V1 = x[7]    # volume in reactor (L)
    # Xv2 = x[8]    # concentration of viable cells in separator(cell/L)
    # Xt2 = x[9]    # total concentration of cells in separator(cell/L)
    # GLC2 = x[10]  # glucose concentration in separator(mM)
    # GLN2 = x[11]  # glutamine concentration in separator(mM)
    # LAC2 = x[12]  # lactate concentration in separator(mM)
    # AMM2 = x[13]  # ammonia concentration in separator(mM)
    # mAb2 = x[14]  # mAb concentration in separator(mg/L)
   ## V2 = x[15]    # volume in separator (L)
    # # temperature as a state -- Ben
    # T = x[16]     #temperature of bioreactor mixture  (d C)
    # # buffer tank
   ## h = x[17]  # Liquid level in the tank (dm)
    # c = x[18]  # Concentration of mAb in the tank (mg/L)
    
    #new state indexing
    # Xv1 = x[0]   # concentration of viable cells in reactor(cell/L)
    # Xt1 = x[1]   # total concentration of cells in reactor(cell/L)
    # GLC1 = x[2]  # glucose concentration in reactor(mM)
    # GLN1 = x[3]  # glutamine concentration in reactor(mM)
    # LAC1 = x[4]  # lactate concentration in reactor(mM)
    # AMM1 = x[5]  # ammonia concentration in reactor(mM)
    # mAb1 = x[6]  # mAb concentration in reactor(mg/L)
    # Xv2 = x[7]    # concentration of viable cells in separator(cell/L)
    # Xt2 = x[8]    # total concentration of cells in separator(cell/L)
    # GLC2 = x[9]  # glucose concentration in separator(mM)
    # GLN2 = x[10]  # glutamine concentration in separator(mM)
    # LAC2 = x[11]  # lactate concentration in separator(mM)
    # AMM2 = x[12]  # ammonia concentration in separator(mM)
    # mAb2 = x[13]  # mAb concentration in separator(mg/L)
    # # temperature as a state -- Ben
    # T = x[14]     #temperature of bioreactor mixture  (d C)
    # # buffer tank
    # c = x[15]  # Concentration of mAb in the tank (mg/L)


    #defining new variables for the normalised states
    x1 = x[0]*(width*xs[0])+lb*xs[0];      x2 = x[1]*(width*xs[1])+lb*xs[1];      x3 = x[2]*(width*xs[2])+lb*xs[2]; 
    x4 = x[3]*(width*xs[3])+lb*xs[3];      x5 = x[4]*(width*xs[4])+lb*xs[4];      x6 = x[5]*(width*xs[5])+lb*xs[5]; 
    x7 = x[6]*(width*xs[6])+lb*xs[6];      x8 = x[7]*(width*xs[7])+lb*xs[7];      x9 = x[8]*(width*xs[8])+lb*xs[8];
    x10 = x[9]*(width*xs[9])+lb*xs[9];     x11= x[10]*(width*xs[10])+lb*xs[10];   x12 = x[11]*(width*xs[11])+lb*xs[11];
    x13 = x[12]*(width*xs[12])+lb*xs[12];  x14= x[13]*(width*xs[13])+lb*xs[13];   x15 = x[14]*(width*xs[14])+lb*xs[14];
    x16 = x[15]*(width*xs[15])+lb*xs[15];  
    
    
    
    #inputs remained the same as defined above
    # x8 = xs[7]; 
    # x16 = xs[15];
    # x18 = xs[17];
    
    # Inputs of upstream model:
    F_in = u[0]    # inlet flow rate (L/min)
    F_1 = u[1]     # outlet flow rate of bioreactor (L/min)
    F_r = u[2]     # u[1]-u[3] # u[2]      # recycle flow rate (L/min)
    F_2 = u[3]     # Outlet flow rate of separator (L/min)
    GLC_in = u[4]  # inlet glucose concentration (mM)
    GLN_in = u[5]  # inlet glutamine concentration (mM)
    Tc = u[6]      # coolant temperature (d C)
    # Inputs of buffer tank:
    F_out_bf = u[7]  # Outlet flow rate of the buffer tank (L/min)  
    F_in_bf = u[3]   # Inlet flow rate of the buffer tank (L/min)
    #c_in_bf = x[14] # Inlet concentration of mAb (mg/L)
    c_in_bf = x14    ##changed to the scaled value of x[14]
    

    #parameters - one of the parameters is state dependent so the parameters had to be brought below the normalised variables
    K_damm = 1.76     # ammonia constant for cell death (mM)
    K_dgln = 9.6e-03/60  # constant for glutamine degradation (min^(-1))
    K_glc = 0.75      # Monod constant for glucose (mM)
    K_gln = 0.075     # Monod constant for glutamine (mM)
    KI_amm = 28.48    # Monod constant for ammonia (mM)
    KI_lac = 171.76   # Monod constant for lactate (mM)
    m_glc = 4.9e-14/60   # maintenance coefficients of glucose (mmol/cell/min)
    n = 2             # constant represents the increase of specific death as ammonia concentration increases (-)
    Y_ammgln = 0.45   # yield of ammonia from glutamine (mmol/mmol)
    Y_lacglc = 2.0    # yield of lactate from glucose (mmol/mmol)
    Y_Xglc = 2.6e8    # yield of cells on glucose (cell/mmol)
    Y_Xgln = 8e8      # yield of cells on glutamine (cell/mmol)
    alpha1 = 3.4e-13/60  # constants of glutamine maintenance coefficient (mM L/cell/min)
    alpha2 = 4        # constants of glutamine maintenance coefficient (mM)
    # mu_max = (0.0016*x[16] - 0.0308)*5/60 # 0.058 * 5  # maximum specific growth rate (min^(-1))
    mu_max = (0.0016*x15 - 0.0308)*5/60 # 0.058 * 5  ##x[16] changed to the scaled value of x16
    #mu_dmax = (-0.0045*x[16] + 0.1682)/60 # 0.06    # maximum specific death rate (min^(-1))
    mu_dmax =0.06/60
    m_mabx = 6.59e-10/60 # constant production of mAb by viable cells (mg/Cell/min) = Qmab_max
    
    # Add new parameters -- Ben
    Delta_H = 5e5     # 4.4e4# 1e-8
    rho = 1560.0      # density of mixture is assumed to be density of glucose [g/L]
    cp = 1.244        # specific heat capacity of mixture is assumed to be that of glucose [J/g/dC]
    U = 4e2
    T_in = 37.0
    
    # Parameters of buffer tank
    D = 5  # Diameter of the tank (dm)
    L = D  # Length of the tank (dm)
    Ac = np.pi*(D/2)**2  # Cross-sectional area (dm^2)
    V = Ac*L  # Volume of the tank (L)
    
    
    # Assumptions: 92% cell recycle rate and 20% mab retention
    Xvr = 0.92 * x1 * F_1 / F_r  # concentration of viable cells in recycle stream(cell/L)
    Xtr = 0.92 * x2 * F_1 / F_r  # total concentration of cells in recycle stream(cell/L)
    GLCr = 0.2 * x3 * F_1 / F_r  # glucose concentration in recycle stream(mM)
    GLNr = 0.2 * x4 * F_1 / F_r  # glutamine concentration in recycle stream(mM)
    LACr = 0.2 * x5 * F_1 / F_r  # lactate concentration in recycle stream(mM)
    AMMr = 0.2 * x6 * F_1 / F_r  # ammonia concentration in recycle stream(mM)
    mAbr = 0.2 * x7 * F_1 / F_r  # mAb concentration in recycle stream (mg/L)
    
    
    # ODEs
    f_lim = (x3 / (K_glc + x3)) * (x4 / (K_gln + x4))
    f_inh = (KI_lac / (KI_lac + x5)) * (KI_amm / (KI_amm + x6))
    mu = mu_max * f_lim * f_inh
    mu_d = mu_dmax / (1 + (K_damm / x6) ** n)
    
    dxdt[0]= (mu * x1 - mu_d * x1 - (F_in / 3.00039503e+03) * x1 + (F_r / 3.00039503e+03) * (Xvr - x1) +w[0]) /(width*xs[0])
    dxdt[1] = (mu * x1 + (F_r / 3.00039503e+03) * (Xtr - x2) - (F_in / 3.00039503e+03) * x2 +w[1]) /(width*xs[1])
    
    Q_glc = mu / Y_Xglc + m_glc
    dxdt[2] = (-Q_glc * x1 + (F_in / 3.00039503e+03) * (GLC_in - x3) + (F_r / 3.00039503e+03) * (GLCr - x3) +w[2]) /(width*xs[2])
    
    m_gln = (alpha1 * x4) / (alpha2 + x4)
    Q_gln = mu / Y_Xgln + m_gln
    dxdt[3] = (-Q_gln * x1 - K_dgln * x4 + (F_in / 3.00039503e+03) * (GLN_in - x4) + (F_r / 3.00039503e+03) * (GLNr - x4) +w[3]) /(width*xs[3])
    
    Q_lac = Y_lacglc * Q_glc
    dxdt[4] = (Q_lac * x1 - (F_in / 3.00039503e+03) * x5 + (F_r / 3.00039503e+03) * (LACr - x5) +w[4]) /(width*xs[4])
    
    Q_amm = Y_ammgln * Q_gln
    dxdt[5] = (Q_amm * x1 + K_dgln * x4 - (F_in / 3.00039503e+03) * x6 + (F_r / 3.00039503e+03) * (AMMr - x6) +w[5]) /(width*xs[5]) 
    dxdt[6] = (m_mabx * x1 - (F_in / 3.00039503e+03) * x7 + (F_r / 3.00039503e+03) * (mAbr - x7) +w[6]) /(width*xs[6])
    
    # dxdt[7] = (F_in + F_r - F_1)/(width*xs[7])
    # dxdt[7] = 0
    
    dxdt[7] = ((F_1 / 2.97694060e+03) * (x1 - x8) - (F_r / 2.97694060e+03) * (Xvr - x8) +w[7]) /(width*xs[7])
    dxdt[8] = ((F_1 / 2.97694060e+03) * (x2 - x9) - (F_r / 2.97694060e+03) * (Xtr - x9) +w[8]) /(width*xs[8])
    dxdt[9] = ((F_1 / 2.97694060e+03) * (x3 - x10) - (F_r / 2.97694060e+03) * (GLCr - x10) +w[9]) /(width*xs[9])
    dxdt[10] = ((F_1 / 2.97694060e+03) * (x4 - x11) - (F_r / 2.97694060e+03) * (GLNr - x11) +w[10]) /(width*xs[10])
    dxdt[11] = ((F_1 / 2.97694060e+03) * (x5 - x12) - (F_r / 2.97694060e+03) * (LACr - x12) +w[11]) /(width*xs[11])
    dxdt[12] = ((F_1 / 2.97694060e+03) * (x6 - x13) - (F_r / 2.97694060e+03) * (AMMr - x13) +w[12]) /(width*xs[12])
    dxdt[13] = ((F_1 / 2.97694060e+03) * (x7 - x14) - (F_r / 2.97694060e+03) * (mAbr - x14) +w[13]) /(width*xs[13])
    # dxdt[15] = (F_1 - F_2 - F_r) /(width*xs[15])
    # dxdt[15] = 0
    
    # Temperature dynamics -- Ben
    # Average mass of an mAb cell is 150 kDa = 2.4908084e-19 g
    dxdt[14] = ((F_in / 3.00039503e+03) * (T_in - x15) + (Delta_H / (rho * cp)) * mu * x1 * 2.4908084e-19 + (U / (3.00039503e+03 * rho * cp)) * (Tc - x15) +w[14]) /(width*xs[14]) 
    
    
    # ODEs of buffer tank:
    dhdt = 1/Ac*(F_in_bf-F_out_bf)
    dcdt = F_in_bf/(Ac*1.50000000e+01)*(c_in_bf-x16)
    # dxdt[17] = dhdt /(width*xs[17])
    # dxdt[17] = 0
    dxdt[15] = (dcdt +w[15]) /(width*xs[15])
                                                                                                                                                         
    #dxdt = [dxdt1, dxdt2, dxdt3, dxdt4, dxdt5, dxdt6, dxdt7, dxdt8, dxdt9, dxdt10, dxdt11, dxdt12, dxdt13, dxdt14, dxdt15, dxdt16, dxdt17, dxdt18]
    return dxdt #(returns scaled values since this is now the scaled version of the ode )


num_measurements=7
num_states=16
def measurement_xy_model(x):            
    C_matrix = np.zeros((num_measurements, num_states))
    C_matrix[0,7]=1
    C_matrix[1,8]=1
    C_matrix[2,9]=1
    C_matrix[3,10]=1
    C_matrix[4,11]=1
    C_matrix[5,12]=1
    C_matrix[6,15]=1
    
    y = np.dot(C_matrix, x)
    return y









